#------------------------------------------------------------------------------
# Import libraries
#------------------------------------------------------------------------------

rm(list=ls())
library(LalRUtils)
libreq(tidyverse, data.table, zoo, tictoc, 
       fst, fixest, PanelMatch, patchwork,
       rio, magrittr, janitor, did, panelView, 
       ggiplot, tictoc, binsreg, interflex)
set.seed(42)
theme_set(lal_plot_theme())

#------------------------------------------------------------------------------





#------------------------------------------------------------------------------
# Define Paths
#------------------------------------------------------------------------------

# R studio
setwd( dirname(rstudioapi::getActiveDocumentContext()$path) )
# R default : unccoment if you use default R
# setwd(getSrcDirectory(function(){})[1])

#------------------------------------------------------------------------------




#------------------------------------------------------------------------------
# Load Data
#------------------------------------------------------------------------------

vcf <- fread("vcf_data_complete.csv", sep = ",")
setnames( vcf, "d", "D")
vcf_data <- copy( vcf )
gfc <- fread("gfc_dta.csv" )


#------------------------------------------------------------------------------



#------------------------------------------------------------------------------
# Plot
#------------------------------------------------------------------------------

# %%
vcf = vcf_data[year>= 1990]
fifth_sched_pre_2k = c("Andhra Pradesh", "Chhattisgarh", 
                       "Gujarat", "Himachal Pradesh",
                       "Orissa", "Rajasthan", "Madhya Pradesh")
pri = vcf[year<=1999 & state %in% fifth_sched_pre_2k]
pri[, first_panch_elec := fcase(
  state == "Andhra Pradesh", 1995,
  state == "Chhattisgarh", 1995,
  state == "Gujarat", 1995,
  state == "Madhya Pradesh", 1994,
  state == "Himachal Pradesh", 1995,
  state == "Orissa", 1997,
  state == "Rajasthan", 1995
)]

# %% flip treatment status since PRI treatment is PESA control and vice vers
pri[, time_pri := year - first_panch_elec]
pri[, D := (1 - sch) * (year >= first_panch_elec)]
pri[, nsch := (1 - sch)]


# %%
pri[, never_treated := max(D) == 0, cellid]
pri[, pref_bin := ntile(cover_1990, 10)]
fitter = function(cut) {
  m = feols(forest_index ~ D | cellid + styear, 
            data = pri[pref_bin >= cut], 
            cluster = ~blk)
  tidy(m)[, 2:3]
}

cutmods = map_dfr(1:10, fitter) %>% setDT
cutmods$n = 1:10
colnames(cutmods)[1:2] = c('beta', 'se')
# %%
(rob_fit_pri = ggplot(cutmods, aes(n, beta)) +
    geom_point(size = 2) +
    geom_errorbar(aes(ymax = beta + 1.96 * se,
                      ymin = beta - 1.96 * se), alpha = 1, width = 0) +
    geom_hline(yintercept = 0) +
    scale_x_continuous(breaks = 1:10) +
    labs(y = "Effect (Forest Index)", 
         x = "Inclusion Threshold Decile of \nForest Cover in 1990",
         caption = '') + 
    theme(legend.position = 'none', 
          axis.text = element_text( size = 22 ), 
          text = element_text( size = 22 ), 
          plot.title = element_text(size = 24))
)
# %%
ggsave( "main_figure7_panelA.pdf", rob_fit_pri, 
        width = 9, height = 5,
        device = cairo_pdf)




# %% one-shot treatment in 2008
vcf[, pref_bin := ntile(cover_1990, 10)]
vcf[, time_fra := year - 2008]
vcf[, D_fra := sch * (year >= 2008)]
fitter = function(cut) {
  m = feols(forest_index ~ D_fra | cellid + styear, 
            data = vcf[pref_bin >= cut], cluster = ~blk)
  tidy(m)[, 2:3]
}

cutmods = map_dfr(1:10, fitter) %>% setDT
cutmods$n = 1:10
colnames(cutmods)[1:2] = c('beta', 'se')
# %%
(FRA_fit_vcf = ggplot(cutmods, aes(n, beta)) +
    geom_point(size = 2) +
    geom_errorbar(aes(ymax = beta + 1.96 * se,
                      ymin = beta - 1.96 * se), 
                  alpha = 1, width = 0) +
    geom_hline(yintercept = 0) +
    scale_x_continuous(breaks = 1:10) +
    labs(y = "Effect (Forest Index)", 
         x= "Inclusion Threshold Decile of \nForest Cover in 1990"
    ) + 
    theme(legend.position = 'none',
          axis.text = element_text( size = 22 ), 
          text = element_text( size = 22 ), 
          plot.title = element_text(size = 24))
  
)
# %%
ggsave( "main_figure7_panelB.pdf", 
        FRA_fit_vcf,
        width = 9, height = 5, 
        device = cairo_pdf)


#------------------------------------------------------------------------------